//system
#include <fstream>
#include <sstream>
#include <vector>
#include <stdexcept>
#include <deque>
#include <string>
#include <iostream>


//boost
#include "boost/tuple/tuple.hpp"

// my headers
#define DEBUG_LEVEL 6
#include <debug_func.h>
#include <blast_results.h>
#include <index_string.h>
#include <portable_timer.h>
#include <codeml_results.h>
#include <progress_indicator.h>
#define _MY_STL_MACROS
#include <stlhelper.h>
#include <ultoa.h>
#include <get_filename.h>
#include <print_error.h>
#include <seq_pairs.h>
#include <open_fstream.h>



#include "get_arg.h"


using std::vector;
using std::pair;
using std::deque;
using std::ofstream;
using std::ostream;
using std::ifstream;
using std::istream;
using std::istringstream;
using std::string;
using std::runtime_error;
using std::cerr;
using std::cin;
using std::cout;
using boost::tie;



struct t_threshold_counts
{
	string					name;
	unsigned				input;
	unsigned            	len_threshold;
	unsigned            	fragment_coverage_threshold;
	unsigned            	residues_coverage_threshold;
	unsigned            	output;
	t_program_args&			args;
	t_threshold_counts(const string& name_, t_program_args& args_):
							name(name_),
							input(0),
							len_threshold(0),
							fragment_coverage_threshold(0),
							residues_coverage_threshold(0),
							output(0),
							args(args_)
							{}
};

enum E_THRESHOLD_GENE_COUNTS{	eINPUT	= 1,
								eLEN	= 2,
								eFRAG	= 4,
								eRES	= 8,
								eOUTPUT	= 16
							};


string left_pad(unsigned num, unsigned len = 8)
{
	string res = ultoa(num);
	if (res.length() < len)
		res = string(len - res.length(), ' ')  + res;
	return res;
}
ostream& operator << (ostream& os, const t_threshold_counts& counts)
{
	os << "\t"		<< left_pad(counts.input)
								<< " " << counts.name <<  " parsed.\n";
	if (counts.args.len_threshold != std::numeric_limits<unsigned>::max())
		os << "\t"	<< left_pad(counts.len_threshold)
								<< " with aligned length >= "
								<< counts.args.len_threshold << "\n";
	if (counts.args.fragment_coverage > 0.0)
		os << "\t" 	<< left_pad(counts.fragment_coverage_threshold)
								<< " with fraction of residues aligned (incl. fragments) > "
								<< counts.args.fragment_coverage << "\n";
	if (counts.args.residues_coverage > 0.0)
		os << "\t" 	<< left_pad(counts.residues_coverage_threshold)
								<< " with fraction of residues aligned (sans fragments) > "
								<< counts.args.residues_coverage << "\n";

	os	 << "\t" 	<< left_pad(counts.output)
								<< " " << counts.name <<  " output.\n";

	return os;
}

//________________________________________________________________________________________

//		count_pairs_over_threshold

//________________________________________________________________________________________
void count_pairs_over_threshold(t_threshold_counts&				cnts,
								const std::deque<unsigned>&		gene_cnts)
{
	for (unsigned i = 0; i < gene_cnts.size(); ++i)
	{
		cnts.input							+= gene_cnts[i] & eINPUT;
		cnts.len_threshold					+= gene_cnts[i] & eLEN;
		cnts.fragment_coverage_threshold	+= gene_cnts[i] & eFRAG;
		cnts.residues_coverage_threshold	+= gene_cnts[i] & eRES;
		cnts.output							+= gene_cnts[i] & eOUTPUT;
	}
	cnts.input							/= eINPUT;
	cnts.len_threshold					/= eLEN;
	cnts.fragment_coverage_threshold	/= eFRAG;
	cnts.residues_coverage_threshold	/= eRES;
	cnts.output							/= eOUTPUT;
}



//________________________________________________________________________________________

//		index_gene_pair

//________________________________________________________________________________________
pair<unsigned, unsigned> index_gene_pair(const string&		name1, 
										 const string&		name2, 
										 t_index_string& 	genes,
										 deque<unsigned>&	gene_cnts)
{
	bool added;
	unsigned gene_index1 = genes.index (name1, added);
	if (added)
		gene_cnts.push_back(1);
	else
		gene_cnts[gene_index1] |= eINPUT;
	unsigned gene_index2 = genes.index (name2, added);
	if (added)
		gene_cnts.push_back(1);
	else
		gene_cnts[gene_index2] |= eINPUT;
	return pair<unsigned,  unsigned>(gene_index1, gene_index2);
}



//________________________________________________________________________________________

//		do_parse_pairs

//________________________________________________________________________________________
struct do_parse_pairs
{
	private:
		string			old_name1;
		string			old_name2;
		t_index_string	genes;
		const unsigned	coverage_required;
	public:
		deque<unsigned>					gene_flags;
		t_threshold_counts&				pair_cnts;
		t_program_args&					args;
		std::string						comment_letters;
		do_parse_pairs(	t_threshold_counts&					pair_cnts_,
						t_program_args&						args_,
						const std::string&					comment_letters_):
				coverage_required(args_.coverage_both ?
												2 :					// both in pair
												1 ),				// either one of pair
				pair_cnts(pair_cnts_),
				args(args_),
				comment_letters(comment_letters_){	}
	void print_error(unsigned cnt_line, const string& prot_name) const
	{
		std_print_error(args.CERR(), "The sequence for \"" + prot_name +
								"\" on line #" + ultoa(cnt_line) +
								" could not be found");
	}
	void operator()(istream& is, unsigned)
	{
		// get data for one protein pair
		string name1, name2;

		// comment if false
		string comment_line;
		if (!seq_pairs::parse_ids(is, name1, name2, comment_line, comment_letters,
								  seq_pairs::ePARSE_IDS_ONLY))
		{						  
			cout	<< comment_line;
			return;
		}		


		bcpp_align::t_align_results aligndata;
		is >> aligndata;

		if (!is)
			throw std::runtime_error("The alignment data is not in the correct format");

		
		// prevent simple repeats
		if (name1 == old_name1 && name2 == old_name2)
			return;
		old_name1 = name1;
		old_name2 = name2;

		// index genes
		unsigned gi1, gi2;
		tie(gi1, gi2) = index_gene_pair(name1, name2, genes, gene_flags);
		++pair_cnts.input;


		// len threshold
		if (aligndata.align_len < args.len_threshold)
			return;
		{
			gene_flags[gi1] |= eLEN;
			gene_flags[gi2] |= eLEN;
			++pair_cnts.len_threshold;
		}



		// fragments threshold
		double min_len = std::min(aligndata.len1, aligndata.len2);
		double max_len = std::max(aligndata.len1, aligndata.len2);
		unsigned pass_threshold = 0;
		if (aligndata.aligned_residues >= args.fragment_coverage * min_len)
		{
			gene_flags[gi1] |= eFRAG;
			gene_flags[gi2] |= eFRAG;
			++pair_cnts.fragment_coverage_threshold;
			++pass_threshold;
		}

		// residues threshold
		if (aligndata.aligned_residues >= args.residues_coverage * max_len)
		{
			gene_flags[gi1] |= eRES;
			gene_flags[gi2] |= eRES;
			++pair_cnts.residues_coverage_threshold;
			++pass_threshold;
		}

		// check if pass coverage threshold (and / or)
		if (pass_threshold < coverage_required )
			return;
		{
			++pair_cnts.output;
			gene_flags[gi1] |= eOUTPUT;
			gene_flags[gi2] |= eOUTPUT;
		}


		cout	<< name1		<< "\t"
				<< name2		<< "\t"
				<< aligndata	<< "\n";

		
	}
};





//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//		Main

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
int main(int argc, char* argv[])
{
	t_program_args args;

	try
	{

	debug_func_cerr(5);

	// don't use C streams
	std::ios::sync_with_stdio(false);


	// get command line arguments
	if (!process_command_line_options(argc, argv, args))
		return 1;




	//
	// filter alignment data
	//
	args.VERBOSE() << "\tReading sequence pair data from [" + 
						(args.file_input.length() ? 
						 get_filename(args.file_input) : 
						 "STDIN") + "]\n";
	t_threshold_counts pair_cnts("Alignment pairs", args);
	do_parse_pairs parse_pairs(pair_cnts, args, args.comment_letters);
	seq_pairs::read_data(parse_pairs, args.file_input, args.VERBOSE());

	// count genes
	t_threshold_counts gene_cnts("identifiers", args);
	count_pairs_over_threshold(gene_cnts, parse_pairs.gene_flags);

	// write summary
	args.VERBOSE() << pair_cnts << "\n" << "\n";
	args.VERBOSE() << gene_cnts << "\n" << "\n";

	// write to file
	if (!args.results_summary.empty())
	{
		ofstream ofs(args.results_summary.c_str());
		if (!ofs)
			throw runtime_error("Could not open " + args.results_summary);
		ofs << pair_cnts << "\n" << "\n";
		ofs << gene_cnts << "\n" << "\n";
	}

	}
	catch (std::exception& e)
	{
		std_print_error(args.CERR(), e.what());
		return 1;
	}

	return 0;
}

